# include <cstdlib>
# include <iostream>
# include <iomanip>
# include <fstream>
# include <cmath>
# include <ctime>

using namespace std;

# include "mpi.h"

int main ( int argc, char *argv[] );
int congruence ( int a, int b, int c, int *error );
int i4_gcd ( int i, int j );
int i4_max ( int i1, int i2 );
int i4_min ( int i1, int i2 );
int i4_sign ( int i );
void lcrg_anbn ( int a, int b, int c, int n, int *an, int *bn );
int lcrg_evaluate ( int a, int b, int c, int x );
int power_mod ( int a, int n, int m );
void timestamp ( );


int main ( int argc, char *argv[] )

{
  int a;
  int an;
  int b;
  int bn;
  int c;
  int error;
  int id;
  int j;
  int k;
  int k_hi;
  int p;
  int u;
  int v;
//
//  Initialize MPI.
//
  MPI::Init ( argc, argv );
//
//  Get the number of processors.
//
  p = MPI::COMM_WORLD.Get_size ( );
//
//  Get the rank of this processor.
//
  id = MPI::COMM_WORLD.Get_rank ( );
//
//  Print a message.
//
  if ( id == 0 )
  {
    timestamp ( );
    cout << "\n";
    cout << "RANDOM_MPI - Master process:\n";
    cout << "  C++ version\n";
    cout << "  The number of processors is P = " << p << "\n";
    cout << "\n";
    cout << "  This program shows how a stream of random numbers\n";
    cout << "  can be computed 'in parallel' in an MPI program.\n";
    cout << "\n";
    cout << "  We assume we are using a linear congruential\n";
    cout << "  random number generator or LCRG, which takes\n";
    cout << "  an integer input and returns a new integer output:\n";
    cout << "\n";
    cout << "    U = ( A * V + B ) mod C\n";
    cout << "\n";
    cout << "  We assume that we want the MPI program to produce\n";
    cout << "  the same sequence of random values as a sequential\n";
    cout << "  program would - but we want each processor to compute\n";
    cout << "  one part of that sequence.\n";
    cout << "\n";
    cout << "  We do this by computing a new LCRG which can compute\n";
    cout << "  every P'th entry of the original one.\n";
    cout << "\n";
    cout << "  Our LCRG works with integers, but it is easy to\n";
    cout << "  turn each integer into a real number between [0,1].\n";
  }

//  A, B and C define the linear congruential random number generator.

  a = 16807;
  b = 0;
  c = 2147483647;

  if ( id == 0 )
  {
    cout << "\n";
    cout << "  LCRG parameters:\n";
    cout << "\n";
    cout << "  A  = " << a << "\n";
    cout << "  B  = " << b << "\n";
    cout << "  C  = " << c << "\n";
  }

  k_hi = p * 10;
//
//  Processor 0 generates 10 * P random values.
//
  if ( id == 0 )
  {
    cout << "\n";
    cout << "  Let processor 0 generate the entire random number sequence.\n";
    cout << "\n";
    cout << "     K    ID         Input        Output\n";
    cout << "\n";

    k = 0;
    v = 12345;
    cout << "  " << setw(4) << k
         << "  " << setw(4) << id
         << "  " << "            "
         << "  " << setw(12) << v << "\n";

    for ( k = 1; k <= k_hi; k++ )
    {
      u = v;
      v = lcrg_evaluate ( a, b, c, u );
      cout << "  " << setw(4) << k
           << "  " << setw(4) << id
           << "  " << setw(12) << u
           << "  " << setw(12) << v << "\n";
    }
  }
//
//  Processor P now participates by computing the P-th part of the sequence.
//
  lcrg_anbn ( a, b, c, p, &an, &bn );

  if ( id == 0 )
  {
    cout << "\n";
    cout << "  LCRG parameters for P processors:\n";
    cout << "\n";
    cout << "  AN = " << an << "\n";
    cout << "  BN = " << bn << "\n";
    cout << "  C  = " << c << "\n";
    cout << "\n";
    cout << "  Have ALL the processors participate in computing\n";
    cout << "  the same random number sequence.\n";
    cout << "\n";
    cout << "     K    ID         Input        Output\n";
    cout << "\n";
  }

  v = 12345;
  for ( j = 1; j <= id; j++ )
  {
    u = v;
    v = lcrg_evaluate ( a, b, c, u );
  }
  k = id;

  cout << "  " << setw(4) << k
       << "  " << setw(4) << id
       << "  " << "            "
       << "  " << setw(12) << v << "\n";
//
//  Now use the "skipping" LCRG to compute the values with indices
//  ID, ID+P, ID+2P, ...,
//
  for ( k = id + p; k <= k_hi; k = k + p )
  {
    u = v;
    v = lcrg_evaluate ( an, bn, c, u );
    cout << "  " << setw(4) << k
         << "  " << setw(4) << id
         << "  " << setw(12) << u
         << "  " << setw(12) << v << "\n";
  }
//
//  Terminate MPI.
//
  MPI::Finalize ( );
//
//  Terminate.
//
  if ( id == 0 )
  {
    cout << "\n";
    cout << "RANDOM_MPI:\n";
    cout << "  Normal end of execution.\n";
    cout << "\n";
    timestamp ( );
  }
  return 0;
}


int congruence ( int a, int b, int c, bool *error )


# define N_MAX 100

  int a_copy;
  int a_mag;
  int a_sign;
  int b_copy;
  int b_mag;
  int b_sign;
  int c_copy;
  int g;
  int k;
  int n;
  float norm_new;
  float norm_old;
  int q[N_MAX];
  bool swap;
  int temp;
  int x;
  int xnew;
  int y;
  int ynew;
  int z;
//
//  Defaults for output parameters.
//
  *error = false;
  x = 0;
  y = 0;
//
//  Special cases.
//
  if ( a == 0 && b == 0 && c == 0 )
  {
    x = 0;
    return x;
  }
  else if ( a == 0 && b == 0 && c != 0 )
  {
    *error = true;
    x = 0;
    return x;
  }
  else if ( a == 0 && b != 0 && c == 0 )
  {
    x = 0;
    return x;
  }
  else if ( a == 0 && b != 0 && c != 0 )
  {
    x = 0;
    if ( ( c % b ) != 0 )
    {
      *error = true;
    }
    return x;
  }
  else if ( a != 0 && b == 0 && c == 0 )
  {
    x = 0;
    return x;
  }
  else if ( a != 0 && b == 0 && c != 0 )
  {
    x = c / a;
    if ( ( c % a ) != 0 )
    {
      *error = true;
    }
    return x;
  }
  else if ( a != 0 && b != 0 && c == 0 )
  {
//  g = i4_gcd ( a, b );
//  x = b / g;
    x = 0;
    return x;
  }
//
//  Now handle the "general" case: A, B and C are nonzero.
//
//  Step 1: Compute the GCD of A and B, which must also divide C.
//
  g = i4_gcd ( a, b );

  if ( ( c % g ) != 0 )
  {
    *error = true;
    return x;
  }

  a_copy = a / g;
  b_copy = b / g;
  c_copy = c / g;
//
//  Step 2: Split A and B into sign and magnitude.
//
  a_mag = abs ( a_copy );
  a_sign = i4_sign ( a_copy );
  b_mag = abs ( b_copy );
  b_sign = i4_sign ( b_copy );
//
//  Another special case, A_MAG = 1 or B_MAG = 1.
//
  if ( a_mag == 1 )
  {
    x = a_sign * c_copy;
    return x;
  }
  else if ( b_mag == 1 )
  {
    x = 0;
    return x;
  }
//
//  Step 3: Produce the Euclidean remainder sequence.
//
  if ( b_mag <= a_mag )
  {
    swap = false;
    q[0] = a_mag;
    q[1] = b_mag;
  }
  else
  {
    swap = true;
    q[0] = b_mag;
    q[1] = a_mag;
  }

  n = 3;

  for ( ; ; )
  {
    q[n-1] = ( q[n-3] % q[n-2] );

    if ( q[n-1] == 1 )
    {
      break;
    }

    n = n + 1;

    if ( N_MAX < n )
    {
      *error = true;
      cout << "\n";
      cout << "CONGRUENCE - Fatal error!\n";
      cout << "  Exceeded number of iterations.\n";
      exit ( 1 );
    }
  }

  y = 0;
  for ( k = n; 2 <= k; k-- )
  {
    x = y;
    y = ( 1 - x * q[k-2] ) / q[k-1];
  }

  if ( swap )
  {
    z = x;
    x = y;
    y = z;
  }

  x = x * a_sign;
//
//  Step 7: Multiply by C, so that X * A + Y * B = C.
//
  x = x * c_copy;
//
//  Step 8: Now force 0 <= X < B.
//
  x = x % b;
//
//  Step 9: Force positivity.
//
  if ( x < 0 )
  {
    x = x + b;
  }

  return x;
# undef N_MAX
}


int i4_gcd ( int i, int j )


{
  int ip;
  int iq;
  int ir;
//
//  Return immediately if either I or J is zero.
//
  if ( i == 0 )
  {
    return i4_max ( 1, abs ( j ) );
  }
  else if ( j == 0 )
  {
    return i4_max ( 1, abs ( i ) );
  }
//
//  Set IP to the larger of I and J, IQ to the smaller.
//  This way, we can alter IP and IQ as we go.
//
  ip = i4_max ( abs ( i ), abs ( j ) );
  iq = i4_min ( abs ( i ), abs ( j ) );
//
//  Carry out the Euclidean algorithm.
//
  for ( ; ; )
  {
    ir = ip % iq;

    if ( ir == 0 )
    {
      break;
    }

    ip = iq;
    iq = ir;
  }

  return iq;
}


int i4_max ( int i1, int i2 )


{
  int value;

  if ( i2 < i1 )
  {
    value = i1;
  }
  else
  {
    value = i2;
  }
  return value;
}


int i4_min ( int i1, int i2 )


{
  int value;

  if ( i1 < i2 )
  {
    value = i1;
  }
  else
  {
    value = i2;
  }
  return value;
}

int i4_sign ( int i )


{
  int value;

  if ( i < 0 )
  {
    value = -1;
  }
  else
  {
    value = 1;
  }
  return value;
}


void lcrg_anbn ( int a, int b, int c, int n, int *an, int *bn )


{
  int am1;
  int anm1tb;
  bool ierror;

  if ( n < 0 )
  {
    cerr << "\n";
    cerr << "LCRG_ANBN - Fatal error!\n";
    cerr << "  Illegal input value of N = " << n << "\n";
    exit ( 1 );
  }

  if ( c <= 0 )
  {
    cerr << "\n";
    cerr << "LCRG_ANBN - Fatal error!\n";
    cerr << "  Illegal input value of C = " << c << "\n";
    exit ( 1 );
  }

  if ( n == 0 )
  {
    *an = 1;
    *bn = 0;
  }
  else if ( n == 1 )
  {
    *an = a;
    *bn = b;
  }
  else
  {
//
//  Compute A^N.
//
    *an = power_mod ( a, n, c );
//
//  Solve
//    ( a - 1 ) * BN = ( a^N - 1 ) mod B
//  for BN.
//
    am1 = a - 1;
    anm1tb = ( *an - 1 ) * b;

    *bn = congruence ( am1, c, anm1tb, &ierror );

    if ( ierror )
    {
      cerr << "\n";
      cerr << "LCRG_ANBN - Fatal error!\n";
      cerr << "  An error occurred in the CONGRUENCE routine.\n";
      exit ( 1 );
    }
  }

  return;
}


int lcrg_evaluate ( int a, int b, int c, int x )


{
  long long int a8;
  long long int b8;
  long long int c8;
  long long int x8;
  int y;
  long long int y8;

  a8 = ( long long int ) a;
  b8 = ( long long int ) b;
  c8 = ( long long int ) c;
  x8 = ( long long int ) x;

  y8 = ( a8 * x8 + b8 ) % c8;

  y = ( int ) ( y8 );

  if ( y < 0 )
  {
    y = y + c;
  }

  return y;
}

int power_mod ( int a, int n, int m )



{
  long long int a_square2;
  int d;
  long long int m2;
  int x;
  long long int x2;

  if ( a < 0 )
  {
    return -1;
  }

  if ( m <= 0 )
  {
    return -1;
  }

  if ( n < 0 )
  {
    return -1;
  }
//
//  A_SQUARE contains the successive squares of A.
//
  a_square2 = ( long long int ) a;
  m2 = ( long long int ) m;
  x2 = ( long long int ) 1;

  while ( 0 < n )
  {
    d = n % 2;

    if ( d == 1 )
    {
      x2 = ( x2 * a_square2 ) % m2;
    }

    a_square2 = ( a_square2 * a_square2 ) % m2;
    n = ( n - d ) / 2;
  }

  while ( x2 < 0 )
  {
    x2 = x2 + m2;
  }

  x = ( int ) x2;

  return x;
}


void timestamp ( void )


{


  static char time_buffer[TIME_SIZE];
  const struct tm *tm;
  size_t len;
  time_t now;

  now = time ( NULL );
  tm = localtime ( &now );

  len = strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm );

  cout << time_buffer << "\n";

  return;

}
